rm(list=ls())

pacman::p_load(lubridate, tidyverse, haven, gtrendsR, rdrobust, rddtools,
               broom, ggpubr, gridExtra)

setwd('put_your_wd_here')

# this code below will re-estimate RDIT results
# for all twitter data across all shootings
# results are written out to a csv file
# shootingPlots <- function(df, date, shooting, term, time=2){
#   
#   cut = lubridate::ymd(date) 
#   start_date = cut - months(time)
#   end_date = cut + months(time)
# 
#   dat <- df %>%
#     filter(date >= start_date & date <= end_date)
#   if(nrow(dat) == 0 | sum(dat$n_tweets) == 0){
#     out = data.frame(Est = 'Insufficient Data',
#                Coef=NA,
#                lower=NA,
#                upper=NA,
#                shooting=shooting,
#                term=term)
#   } else {
#   dat$running_day <- as.numeric(dat$date - cut)
#   sd_est <- sd(dat$n_tweets, na.rm=T)
#   house_rdd <- rdd_data(y=dat$n_tweets, x=dat$running_day, cutpoint=0)
#   bw_ik <- tryCatch({rdd_bw_ik(house_rdd)}, error = function(e){NA})
#   if(is.nan(bw_ik)){
#     reg_nonpara <- NA
#   } else {
#   reg_nonpara <- tryCatch({rdd_reg_np(rdd_object=house_rdd, bw=bw_ik)}, 
#     error=function(e){NA})
#   }
#   if(!is.na(reg_nonpara[1])){
#   out <- data.frame(Est = 'Rdd Tools Coef',
#                     Coef=reg_nonpara$coefficients/sd_est,
#                     lower=(reg_nonpara$coefficients - 1.96 * reg_nonpara$coefMat[2])/sd_est,
#                     upper=(reg_nonpara$coefficients + 1.96 * reg_nonpara$coefMat[2])/sd_est,
#                     shooting=shooting,
#                     term=term)
#   } else {
#     out = data.frame(Est = 'Insufficient Data',
#                      Coef=NA,
#                      lower=NA,
#                      upper=NA,
#                      shooting=shooting,
#                      term=term)
#     }
#   }
#   return(out) 
# }

# load datasets
# guncontrol <- read_csv('datasets/_full_gun_control.csv')
# gunrights <- read_csv('datasets/_full_gun_rights.csv')
# nra <- read_csv('datasets/_full_nra.csv')
# everytown <- read_csv('datasets/_full_everytown.csv')
# recycling <- read_csv('datasets/_full_recycling.csv')
# 
# everytown$date <- mdy(everytown$date)

# load list of shootings
allshootings <- read_csv('datasets/list_of_shootings.csv',col_select = 1:4)
allshootings$date <- mdy(allshootings$date)
include <- allshootings %>%
  filter(include==1) %>%
  slice(1:19)

# makePlots <- function(dat, term){
#   listSave <- vector('list', nrow(include))
#   for(i in 1:nrow(include)){
#     listSave[[i]] <- tryCatch({shootingPlots(df=dat,
#                   date=include$date[i], 
#                   shooting=include$shooting[i],
#                   term=term)}, error=function(e){NA})
#     print(i)
#   }
#   listSave <- bind_rows(listSave)
#   return(listSave)
# }
# 
# dat <- bind_rows(
#   makePlots(nra, '#NRA'),
#   makePlots(everytown, '#everytown'),
#   makePlots(guncontrol, '#guncontrol'),
#   makePlots(gunrights, '#gunrights'),
#   makePlots(recycling, '#recycling')
# )
# 
# # save results
# write_csv(dat, file='datasets/rdit_results_twitter.csv')

# read in results
dat <- read_csv('datasets/rdit_results_twitter.csv')

# make plots
dat <- dat %>%
  mutate(shooting = factor(shooting, levels=rev(c(
    include$shooting
  )))) %>%
  mutate(term = factor(term, levels=c(
    '#guncontrol','#gunrights','#everytown','#NRA','#abortion','#metoo','#recycling'
  ))) %>%
  na.omit()

merge_expand <- expand.grid(dat$shooting,dat$term) %>%
  unique() %>%
  rename(shooting=Var1, term=Var2)

dat <- left_join(merge_expand, dat, by=c('term','shooting'))

# META ANALYSIS
calcMeta <- function(data, name, outcome, term){
  m.gen <- meta::metagen(TE = Coef,
                   seTE = se,
                   studlab = shooting,
                   data = data,
                   sm = "SMD",
                   fixed = FALSE,
                   random = TRUE)
  results <- tibble(shooting = name, 
                    Coef = m.gen$TE.random,
                    lower = m.gen$lower.random,
                    upper = m.gen$upper.random,
                    outcome = outcome,
                    term=term)
  return(results)
}

dat_meta <- dat %>%
  mutate(se=(upper - Coef)/1.96,
         top10 = ifelse(shooting %in% allshootings$shooting[allshootings$top10==1],1,0))

res_meta <- bind_rows(
  calcMeta(dat_meta[dat_meta$term == '#NRA',],name='Meta Estimate (Top 20)',outcome = '',term = 'Meta #NRA'),
  calcMeta(dat_meta[dat_meta$term == '#everytown',],name='Meta Estimate (Top 20)',outcome = '',term = 'Meta #everytown'),
  calcMeta(dat_meta[dat_meta$term == '#NRA' & dat_meta$top10==1,],name='Meta Estimate (Top 10)',outcome = '',term = 'Meta #NRA'),
  calcMeta(dat_meta[dat_meta$term == '#everytown' & dat_meta$top10==1,],name='Meta Estimate (Top 10)',outcome = '',term = 'Meta #everytown'),
  calcMeta(dat_meta[dat_meta$term == '#guncontrol',],name='Meta Estimate (Top 20)',outcome = '',term = 'Meta #guncontrol'),
  calcMeta(dat_meta[dat_meta$term == '#gunrights',],name='Meta Estimate (Top 20)',outcome = '',term = 'Meta #gunrights'),
  calcMeta(dat_meta[dat_meta$term == '#guncontrol' & dat_meta$top10==1,],name='Meta Estimate (Top 10)',outcome = '',term = 'Meta #guncontrol'),
  calcMeta(dat_meta[dat_meta$term == '#gunrights' & dat_meta$top10==1,],name='Meta Estimate (Top 10)',outcome = '',term = 'Meta #gunrights'),
  calcMeta(dat_meta[dat_meta$term == '#recycling',],name='Meta Estimate (Top 20)',outcome = '',term = 'Meta #recycling'),
  calcMeta(dat_meta[dat_meta$term == '#recycling' & dat_meta$top10==1,],name='Meta Estimate (Top 10)',outcome = '',term = 'Meta #recycling')
)

####
# calulate meta analyses for race analysis
####

race_est <- read_csv('datasets/shooting_race.csv')
race_est <- race_est %>%
  filter(shooting %in% dat_meta$shooting)
terciles <- quantile(race_est$mean.white, na.rm=T, c(0,0.33,0.66,1))
race_est$tercile <- ifelse(race_est$mean.white < terciles[2],1,0)
race_est$tercile <- ifelse(race_est$mean.white > terciles[2] & race_est$mean.white < terciles[3],2,race_est$tercile)
race_est$tercile <- ifelse(race_est$mean.white > terciles[3],3,race_est$tercile)
datnew <- left_join(dat_meta, race_est, by='shooting')

meta_top3 <- calcMeta(datnew %>% mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 3 & term == '#NRA'),  
                      'Meta Estimate (Top)','Meta Analysis Estimate','Meta Analysis')
meta_top2 <- calcMeta(datnew %>% mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 2 & term == '#NRA'), 
                      'Meta Estimate (Middle)','Meta Analysis Estimate','Meta Analysis')
meta_top1 <- calcMeta(datnew %>% mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 1 & term == '#NRA'), 
                      'Meta Estimate (Lower)','Meta Analysis Estimate','Meta Analysis')
meta_tercile_NRA <- bind_rows(meta_top3, meta_top2, meta_top1) %>%
  mutate(category = '#NRA')

meta_top3 <- calcMeta(datnew %>% mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 3 & term == '#everytown'),  
                      'Meta Estimate (Top)','Meta Analysis Estimate','Meta Analysis')
meta_top2 <- calcMeta(datnew %>% mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 2 & term == '#everytown'), 
                      'Meta Estimate (Middle)','Meta Analysis Estimate','Meta Analysis')
meta_top1 <- calcMeta(datnew %>% mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 1 & term == '#everytown'), 
                      'Meta Estimate (Lower)','Meta Analysis Estimate','Meta Analysis')
meta_tercile_everytown <- bind_rows(meta_top3, meta_top2, meta_top1) %>%
  mutate(category = '#everytown')

meta_top3 <- calcMeta(datnew %>% mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 3 & term == '#guncontrol'),  
                      'Meta Estimate (Top)','Meta Analysis Estimate','Meta Analysis')
meta_top2 <- calcMeta(datnew %>% mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 2 & term == '#guncontrol'), 
                      'Meta Estimate (Middle)','Meta Analysis Estimate','Meta Analysis')
meta_top1 <- calcMeta(datnew %>% mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 1 & term == '#guncontrol'), 
                      'Meta Estimate (Lower)','Meta Analysis Estimate','Meta Analysis')
meta_tercile_guncontrol <- bind_rows(meta_top3, meta_top2, meta_top1) %>%
  mutate(category = '#guncontrol')

meta_top3 <- calcMeta(datnew %>% mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 3 & term == '#gunrights'),  
                      'Meta Estimate (Top)','Meta Analysis Estimate','Meta Analysis')
meta_top2 <- calcMeta(datnew %>% mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 2 & term == '#gunrights'), 
                      'Meta Estimate (Middle)','Meta Analysis Estimate','Meta Analysis')
meta_top1 <- calcMeta(datnew %>% mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 1 & term == '#gunrights'), 
                      'Meta Estimate (Lower)','Meta Analysis Estimate','Meta Analysis')
meta_tercile_gunrights <- bind_rows(meta_top3, meta_top2, meta_top1) %>%
  mutate(category = '#gunrights')

meta_res_race <- bind_rows(meta_tercile_NRA, meta_tercile_everytown, meta_tercile_gunrights, meta_tercile_guncontrol)
write_csv(meta_res_race, 'datasets/metaanalysis/meta_race_tweets.csv' )

# res meta lower 20
resmeta2 <- read_csv('datasets/lower_20_tweets.csv')
res_meta <- bind_rows(res_meta, resmeta2)
res_meta <- res_meta %>%
  mutate(term = case_when(
    term == '#NRA' ~ "Meta #NRA",
    term == '#everytown' ~ "Meta #everytown",
    term == '#guncontrol' ~ "Meta #guncontrol",
    term == '#gunrights' ~ "Meta #gunrights",
    term == '#recycling' ~ "Meta #recycling",
    TRUE ~ term
  ))

res_meta %>%
  mutate(z = 'Twitter All') %>%
  write_csv(file='datasets/meta_twitter_all.csv')

dat <- bind_rows(dat, res_meta) %>%
  mutate(statsig = ifelse((lower > 0 & upper > 0) | (lower < 0 & upper < 0),1,0))
dat_ideas <- dat %>% 
  filter(term %in% c('#guncontrol','#gunrights','Meta #gunrights',
                     'Meta #guncontrol'))
dat_gc <- dat %>% 
  filter(term %in% c('#everytown','#NRA','Meta #everytown',
                     'Meta #NRA'))
dat_gr <- dat %>% 
  filter(term %in% c('#recycling','Meta #recycling'))

plot_ideas <- dat_ideas %>%
  mutate(shooting = factor(shooting, levels=rev(c(dat_ideas$shooting[20:38],
                                                  'Meta Estimate (Top 10)',
                                                  'Meta Estimate (Top 20)',
                                                  'Meta Estimate (Bottom 20)'))),
         term = factor(term, levels=rev(c(
           '#guncontrol','#gunrights','Meta #guncontrol','Meta #gunrights')
           ))) %>%
  ggplot(aes(x=shooting, y=Coef, ymin=lower, ymax=upper, color=factor(statsig), shape=term)) +
  geom_errorbar(width=0, position=position_dodge(width=0.65)) +
  geom_point(position=position_dodge(width=0.65), size=2, fill='white') +
  coord_flip() +
  geom_hline(yintercept=0, linetype=2, color='red') +
  labs(title='(A) Policy Positions ', x='',y='') +
  scale_color_manual(values=c('grey50','black'), guide='none') +
  theme_bw() +
  scale_shape_manual(values=c(17,16,24,21), name='') +
  scale_y_continuous(limits=c(-3, 6)) +
  theme(legend.position='none', 
        plot.title = element_text(hjust = 0.5))
plot_ideas

dat_gc <- dat_gc %>% mutate(black = ifelse((lower < 0 & upper < 0) | 
                                             (lower > 0 & upper > 0) ,1,0))
plot_gc <- dat_gc %>%
  mutate(shooting = factor(shooting, levels=rev(c(dat_ideas$shooting[20:38],
                                                  'Meta Estimate (Top 10)',
                                                  'Meta Estimate (Top 20)',
                                                  'Meta Estimate (Bottom 20)'))),
         term = factor(term , levels=c(
           '#NRA','#everytown','Meta #NRA','Meta #everytown'
         ))) %>%
  ggplot(aes(x=shooting, y=Coef, ymin=lower, ymax=upper, color=factor(black), shape=term, group=term)) +
  theme_bw()

plot_gc <- plot_gc + 
  coord_flip()  +
  geom_errorbar(width=0, position=position_dodge(width=0.65)) +
  geom_point(position=position_dodge(width=0.65), size=2, fill='white') +
  geom_hline(yintercept=0, linetype=2, color='red') +
  labs(title='(B) Gun Orgs', x='',y='') +
  scale_color_manual(values=c('grey50','black'), guide='none')  +
  scale_shape_manual(values=c(24, 21, 17, 16), name='') +
  scale_y_continuous(limits=c(-3, 6)) + 
  theme(legend.position='none', 
        plot.title = element_text(hjust = 0.5),
        panel.grid.major = element_line(size = 0.5, linetype = 'solid'),
        axis.text.y=element_blank()) 
plot_gc 

plot_gr <- dat_gr %>%
  mutate(term = factor(term, levels=c(
     '#recycling',
     'Meta #recycling'
  )),
  shooting = factor(shooting, levels=rev(c(dat_gr$shooting[1:19],'Meta Estimate (Top 10)',
                                           'Meta Estimate (Top 20)','Meta Estimate (Bottom 20)')))) %>%
  ggplot(aes(x=shooting, y=Coef, ymin=lower, ymax=upper, color=as.factor(statsig), shape=factor(term))) +
  coord_flip() +
  geom_hline(yintercept=0, linetype=2, color='red') +
  labs(title='(C) Placebo', x='',y='') +
  theme_bw() +
  geom_errorbar(width=0, 
                position=position_dodge(width=0.25)) +
  geom_point(size=2, fill='white')  +
  scale_shape_manual(values=c(21, 16)) +
  scale_color_manual(values=c('grey50','Black'),guide='none') +
  theme(legend.position='none',   
        plot.title = element_text(hjust = 0.5),
        axis.text.y = element_blank()) +
  scale_y_continuous(limits=c(-3, 6))
plot_gr

plots <- ggarrange(plot_ideas, plot_gc,plot_gr,
                   heights = c(7,7,7),
                   widths = c(8,5,5),
                   ncol = 3, nrow = 1)
plots
plots <- annotate_figure(plots, bottom = text_grob("Effect on Tweets\n(Standard Deviations)",family='Helvetica'))
ggsave(plots, width=10, height=5, file='tweets_all.png', bg='white')


# Lower 20 ----
include <- allshootings %>%
  filter(include==0)

dat_new <- bind_rows(
  makePlots(nra, '#NRA'),
  makePlots(everytown, '#everytown'),
  makePlots(guncontrol, '#guncontrol'),
  makePlots(gunrights, '#gunrights'),
  makePlots(recycling, '#recycling')
)

dat_meta <- dat_new %>% mutate(se=(upper - Coef)/1.96)
res_meta <- bind_rows(
  calcMeta(dat_meta[dat_meta$term == '#NRA',],name='Meta Estimate (All)',outcome = '',term = '#NRA'),
  calcMeta(dat_meta[dat_meta$term == '#everytown',],name='Meta Estimate (All)',outcome = '',term = '#everytown'),
  calcMeta(dat_meta[dat_meta$term == '#guncontrol',],name='Meta Estimate (All)',outcome = '',term = '#guncontrol'),
  calcMeta(dat_meta[dat_meta$term == '#gunrights',],name='Meta Estimate (All)',outcome = '',term = '#gunrights'),
  calcMeta(dat_meta[dat_meta$term == '#recycling',],name='Meta Estimate (All)',outcome = '',term = '#recycling'),
)

res_meta2 <- bind_rows(
  calcMeta(dat_meta[dat_meta$term == '#NRA',],name='Meta Estimate (Bottom 20)',outcome = '',term = '#NRA'),
  calcMeta(dat_meta[dat_meta$term == '#everytown',],name='Meta Estimate (Bottom 20)',outcome = '',term = '#everytown'),
  calcMeta(dat_meta[dat_meta$term == '#guncontrol',],name='Meta Estimate (Bottom 20)',outcome = '',term = '#guncontrol'),
  calcMeta(dat_meta[dat_meta$term == '#gunrights',],name='Meta Estimate (Bottom 20)',outcome = '',term = '#gunrights'),
  calcMeta(dat_meta[dat_meta$term == '#recycling',],name='Meta Estimate (Bottom 20)',outcome = '',term = '#recycling'),
)

write_csv(res_meta2, file='datasets/lower_20_tweets.csv')

twit_low_attn <- dat_new %>%
  bind_rows(res_meta) %>%
  mutate(stat_sig = ifelse((lower < 0 & upper < 0) | (lower > 0 & upper > 0),1,0),
         shooting = factor(shooting, levels=rev(c(
           allshootings$shooting, 'Meta Estimate (All)'
         ))),
         term = factor(term , levels=c(
           '#guncontrol','#gunrights','#everytown','#NRA','#recycling'
         ))) %>%
  ggplot(aes(shooting, Coef, ymin=lower, ymax=upper, color=factor(stat_sig))) +
  geom_hline(yintercept=0, linetype=2, color='red') +
  geom_errorbar(width=0.1) +
  geom_point(size=2, shape=21, fill='white') +
  facet_wrap(~term, ncol=5) +
  theme_bw() + 
  labs(x='',y='RDIT Estimate') +
  coord_flip() +
  scale_color_manual(values=c('grey50','black'), guide='none') 
ggsave(twit_low_attn, width=8, height=4, filename = 'figures/twit_low_attn.png')

